AD  750185 


ANALYSIS  OF  JET  EXHAUSTS 
WITH  REPETITIVE 
SHOCK  STRUCTURE 

FINAL  REPORT 
September  1972 


D  D  C 


Reproduced  by 

NATIONAL  TECHNICAL 
INFORMATION  SERVICE 


U  5  Deportment  of  Commerce 
Springfield  VA  22151 


X 


HREC-1251-2 
LMSC-HREC  D306224 


LOCKHEED  MISSILES  A  SPACE  COMPANY  INC. 
HUNTSVILLE  RESEARCH  A  ENGINEERING  CENTER 
HUNTSVILLE  RESEARCH  PARK 
4800  BRADFORD  DRIVE,  HUNTSVILLE,  ALABAMA 


D  D  C 

I“,u  »  m 


ANALYSIS  OF  JET  EXHAUSTS 
WITH  REPETITIVE 
SHOCK  STRUCTURE 

FINAL  REPORT 
September  1972 


Contract  DAAH01 -71 -C - 125 1 


Prepared  for  Commanding  General,  U.  S.  Army  Missile  Command 
Redstone  Arsenal,  Alabama  3  5812 


by 

R. J.  Prozan 
A.  W.  Ratliff 


APPROVED:  'Jri£  ^ 

^  Mf  D _ K-L 


John  W.  Benefield/ Supervisor 
Fluid  Mechanics  Section 


Aeromechanics  Department 


.  1  A 


Approved  for  public  releases 
Distribution  Unlimited  # 


J.S.  Farrior 
Resident  Director 


LMSC-HREC  D  306 224 


FOREWORD 

This  document  reports  the  results  of  an  analysis 
performed  by  Lockheed's  Huntsville  Research  &  Engi¬ 
neering  Center  under  Contract  DAAH01-71-C-1251  for 
the  Army  Missile  Command,  Redstone  Arsenal,  Alabama. 
The  program  was  monitored  by  Mr.  H,  Tracy  Jackson 
of  the  Missile  Systems  Laboratory. 
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SUMMARY 

Described  herein  is  an  advancement  in  the  state  of  the  art  in  calculation 
of  the  mixing  of  two  concentric  (or  parallel)  streams.  The  analysis  treats  the 
equations  of  motion  governing  flow  of  subsonic  or  supersonic  jets  mixing  with 
the  surrounding  atmosphere. 

Several  analyses,  notably  one  of  AeroChem  (Ref.  1)  treat  the  mixing  of 
concentric  streams.  These  analyses,  however,  assume  that  there  is  no  lateral 
pressure  gradient.  This  apriori  assumption  restricts  these  analyses  to  balanced 
jets  which  do  not  vigorously  mix  or  react. 

An  exhaustive  study  of  full  Navier -Stokes  solutions  and/or  simplified 
mathematical  models  which  retain  an  elliptic  nature  did  not  produce  an 
economically  feasible  approach  to  the  solution  of  the  jet  mixing  with  shock 
waves.  The  approach  discussed  here  was  selected  as  an  alternative  since 
inclusion  of  the  radial  momentum  equation  and  the  treatment  of  lateral  pres¬ 
sure  gradients  is  an  improvement  over  current  production  calculations  which 
may  be  achieved  with  nominal  computer  run  times  and  storage. 

The  governing  equations  of  motion  are  discussed  along  with  the  numerical 
analog  and  the  method  of  solution.  Results  of  sample  calculations  are  presented 
and  discussed.  An  appendix  attached  to  the  report  discusses  the  computer  pro¬ 
gram  and  contains  an  input/output  guide. 

Although  the  solution  presented  is  not  the  final  answer  to  the  jet  mixing 
problem  it  is  a  necessary  step  in  the  evolution  of  numerical  solutions  to  this 
complex  flow  system. 
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NOMENCLATURE 


Symbol 

C 

Pi 

EI,  II,  IU,  IV,  V 
f 

G 

g 

H 

h. 

1 

P 

q 

R 

S,  n 
T 

WI,  II,  III,  IV,  V 

y 

Greek 

a 

6 

9 

P 

# 

* 

P 

o 

* 

b) 


Description 

species  specific  heat  at  constant  pressure 

Error  measure  defined  in  text 
any  function 

system  positive  definite  error  measure 

local  positive  definite  error  measure 

total  enthalpy 

species  static  enthalpy 

pressure 

streamwise  velocity 
universal  gas  constant 
streamwise  and  normal  coordinates 
temperature 

weighting  or  scaling  functions  defined  in  text 
distance  from  centerline 

species  mass  fraction 
step  modifier 
flow  angle 

eddy  viscosity  coefficient 
any  independent  variable 
density  of  the  material 

0  for  two-dimensional  flow;  1  for  axisymmetric 
molecular  weight 
species  production  rate 
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Section  1 
INTRODUCTION 

The  conventional  jet  engine  expels  a  stream  of  hot  combustion  products 
into  the  surrounding  atmosphere.  These  gases  radiate  in  the  visible  as  well 
as  the  infrared  spectrum.  The  resultant  radiation  may  be  detected  by  suitable 
instrumentation  and  used  to  track  the  vehicle. 

From  a  military  point  of  view  this  has  an  impact  on  defensive  weapons 
design  (detection  and  aiming)  as  well  as  offensive  weapons  (protection  against 
detection  and  aiming  by  hostile  forces).  The  capability  to  predict  the  gas 
dynamic  behavior  and  salient  features  of  the  jet/atmosphere  interaction  is 
therefore  an  important  part  of  the  overall  problem. 

The  flowfield  structure  is  very  complex  and  has  occupied  the  efforts  of 
many  researchers  over  the  years.  Considerable  effort  has  been  devoted  to 
each  of  several  important  regions  within  the  exhaust  such  as  the  near  field  or 
base  region,  the  shear  layer,  the  far  field  fully  mixed  region  and,  in  the  case 
of  supersonic  jets,  the  internal  shock  structure  of  the  inviscid  region  of  the  jet. 

Mixed  flow  (elliptic,  hyperbolic)  exists  in  these  exhaust  structures  which 
tremendously  complicates  the  analysis.  Successful,  that  is  practical,  analysis 
must  forward  march  from  the  jet  exit  downstream  in  solving  the  equations  of 
motion  and  therefore  are  of  the  parabolic  or  hyperbolic  type.  The  region  of 
interest  is  so  large  and  the  resolution  requirements  so  small  that  an  elliptic 
solution  becomes  impractical. 

In  this  study,  many  methods  were  investigated  in  the  hopes  that  the  prob¬ 
lems  could  be  solved  in  an  economical  manner.  Some  of  the  numerical  tech¬ 
niques  investigated  were  asymptotic  time -dependent  formulations  and  relaxation 
solutions  using  finite  difference  and  finite  element  formulations  for  viscous  and 
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inviscid  formulations.  No  feasible  method  was  found  which  provided  the 
necessary  resolution  over  the  entire  region  of  interest  while  treating  the 
problem  as  elliptic.  Machine  times  must  be  measured  in  hours,  while 
storage  requirements  range  into  hundreds  of  thousands  of  words. 

A  discussion  of  the  state  of  the  art  of  mixing  analysis  "Review  of  Eddy 
Viscosity  Models  for  Jet  Engine  Exhaust/Air  Mixing"  by  B.  J.  Audeh,  JLMSC- 
HREC  D225588,  June  1972  was  also  produced  during  the  course  of  this  study. 
The  methods  discussed  are  all  parabolic  because  the  parabolic  analysis  offers 
the  best  avenue  to  provide  engineering  solutions  to  these  problems.  Most  if 
not  all  such  analyses  assume  a  balanced  jet;  i.e.,  no  lateral  pressure  gradients. 
It  was  possible  to  remove  this  restriction,  however,  and  thereby  strengthen 
the  capability  to  study  jet  exhaust/  atmosphere  interaction  problems.  The 
following  discussion  is  concerned  with  the  technical  approach  to,  and  details 
of,  a  forward  marching  solution  to  the  equations  of  motion  with  lateral  pressure 
gradients . 
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Section  2 

TECHNICAL  DISCUSSION 

The  analysis  of  a  jet  exhausting  into  the  atmosphere  is  truly  a  complex 
procedure.  Certain  calculations  have  been  produced  which  under  conditions 
consistent  with  the  basic  assumptions  have  served  the  designer  and  analyst 
well.  For  supersonic  inviscid  flow  the  method  of  characteristics  may  be  em¬ 
ployed  to  produce  a  rapid,  accurate  description  of  the  resultant  flow  field  until 
a  Mach  disk  is  encountered.  This  type  of  analysis  is  generally  most  useful  for 
highly  underexpanded  jets  in  which  the  Mach  disk  is  many  exit  radii  downstream 
and  far  beyond  the  region  of  interest. 

At  the  other  end  of  the  spectrum,  a  balanced  jet  may  be  reasonably 
handled  by  parabolic  analyses,  a  family  of  solutions  to  which  most  mixing 
programs  belong. 

Perhaps  the  most  complex  condition  is  the  moderately  expanded  jet  in 
which  both  the  lateral  pressure  gradients  and  the  mixing  at  the  periphery  of 
the  jet  are  of  equal  importance.  It  is  technically  feasible  to  employ  one  or 
more  of  the  many  asymptotic  time -dependent  solutions  or  elliptic  relaxation 
procedures  to  analyze  such  fields.  Several  approaches  were  investigated  during 
the  course  of  this  study  leading  to  the  conclusion  that  the  technical  feasibility 
is  severely  compromised,  however,  by  the  economics  of  computer  solutions  of 
this  type.  For  the  problems  of  interest  to  this  study,  one  must  be  willing  to 
accept  several  hours  of  computation  time  and  hundreds  of  thousands  of  words 
of  storage  on  the  most  advanced  machine  in  order  to  perform  such  a  calculation. 

The  approach  taken  here  is  to  treat  the  moderately  expanded  jet  as  a 
forward -marching  (parabolic)  analysis,  yet  retain  the  lateral  momentum 
equation.  This  procedure  yields  a  solution  which  is  rapid,  yet  retains  the 
important  features  of  the  flow  field. 
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2.1 

EQUATIONS  OF  MOTION 

The  analysis  begins  by  considering  the  equations  of  motion  normally 
used  in  balanced  jet  mixing  analyses  but  in  addition  considers  the  lateral 
momentum  equation.  These  equations  presented  in  a  streamwise  normal 
coordinate  system  may  be  found  in  many  texts  and  in  the  interest  of  brevity 
are  merely  restated  here; 

• 

Global  Continuity 

(1) 

• 

Streamwise  Momentum 

pqia  =  .  3E  +  uji^a  +  CT  CO80-  •&! 

P,8S  9S  P(3n2  f  *»l 

(2) 

• 

Normal  Momentum 

2  80  9d  - 

M  8S  +  ta  =  0 

(3) 

• 

Energy 

_  8H  ..  |a2H  .  _  cos 9  8h| 

0“  8S  *  <*(^2  +0  y  9n| 

(4) 

• 

Species  Continuity 

.  tai  ..1  ^“i  .  cosQ  da)  .  • 

<*  as  =  +°  y  anj  +  ui 

(5) 

• 

State 

P  •  pRT  2  T 
i  =  l 

(6) 
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where  I 

h  =  2  «i  hi +  ‘i2/2 
i=l 

and  the  Lewis  and  Prandtl  numbers  have  been  taken  as  unity. 

Thermochemical  information  exists  so  that  for  any  species 

h.  =  h.(T) 

is  known.  It  will  be  convenient,  however,  to  refer  to  the  "specific  heat"  for 
a  species; 

9h. 

cp.  =  “aT  • 

A  more  convenient  form  of  the  energy  equation  is  to  express  the  relations 
explicitly  in  terms  of  temperature.  Recalling  the  definition  of  enthalpy  H,  after 
some  rearrangement,  the  following  can  be  obtained: 


From  the  species  continuity  equation  it  is  evident  that  species  distribution  is 
altered  due  to  the  mixing  and  to  species  production  in  the  event  that  chemical 
reactions  are  of  interest.  An  excellent  finite  rate  computation  produced  by 
AeroChem  (Ref.  2)  has  been  used  in  this  analyses  to  calculate  the  species 
production  term.  The  reader  is  referred  to  Refs.  1  and  2  for  details  of  the 
thermochemical  calculation.  By  setting  u  to  zero  a  frozen  analysis  will  result. 
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After  the  basic  partial  differential  equations  which  have  been  chosen  to 
represent  the  physical  situation  have  been  written,  a  procedure  must  be  de¬ 
veloped  which  solves  the  equations  for  the  flow  parameters.  Since  a  closed 
form  solution  is  unknown,  one  must  resort  to  a  numerical  procedure. 

2.2  NUMERICAL  ANALOG 

A  schematic  of  an  idealized  coaxial  mixing  region  is  shown  in  Fig.  1. 


r*~uc”H 


Initial 

Grid 

Line 


Potential  Core 


Fig.  1  -  Idealized  Coaxial  Mixing 
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An  initial  grid  distribution  is  chosen  and  a  step  is  taken  in  the  axial  direction 
away  from  each  grid  point,  c  reating  a  new  grid  line.  A  typical  section  of  the 
grid  is  shown  in  Fig.  2.  Referring  to  this  figure  the  derivatives  of  a  function 


Fig.  2  -  Typical  Grid  Section 
will  be  approximated  in  the  following  manner: 


8f  ~  ^m+l.n+1  "  *m,n-l 

df  ~  *m+l,n+l 

*m-l,  n+1 

as  _  AS 

911  ym+l,n+I 

"  ym- 1,  n+1 

92f 

^m+l.n+l  +*m-l,n+l  "  9*m,n+l 

9n2  =  ' 

(ym+l,n+l  ‘  ym-l,n+l) 

In  addition,  all  zeroth-order  derivatives  of  f  appearing  in  the  governing  equa 
tions  will  be  approximated  by  the  value  of  the  function  at  m,  n+1. 


As  previously  mentioned,  the  analysis  will  be  parabolic  and  this  means 
that,  while  station  n  is  known  in  its  entirety,  station  n+1  points  must  all  be 
v-’.ved  simultaneously.  In  order  to  accomplish  this  solution  of  a  set  of  non- 
algebraic  equations,  the  error  minimization  technique  (Ref.  3  )  will  be 
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employed.  This  technique  creates  a  positive -definite  function  of  the  errors 
which  exist  in  the  flowfield  equations  and  systematically  reduces  these  errors 
until  an  acceptable  level  (convergence)  is  reached. 

Refer  to  Fig.  2.  For  simplicity,  let  the  index  m,  n  be  represented  to 
zero  and  the  other  points  as  indicated  by  the  appropriate  integers.  The 
numerical  analogs  then  become 


Global  Continuity 


P2q2 


(H)2  +  P2(H)2  +  q2  (H)2  +  °  P2q2 


=  E, 


Streamwise  Momentum 


P2q2 


(»a\  +(9e)  .  |/3fa.\  i  „ coi92 (S3.)  I  . 
M2+W2  11  y2  W/2j 


TI 


(8) 


(9) 


Normal  Momentum 


P2q2 


/ae\  ,  (dp\ 
\*)+\9nL  ~ 


III 


Energy 


C 


2  2  p 


O  /£r\  nr  l$-l\  CQ892  /3T\  ..  /  0a\ 

-  “2  (9s)2  -  “  cp  (8n2)2  +  °  y2  U)2  •"U); 

•twMELto&L  ■  E:v 

•  1  fai  l  *  “ 


f 

2 


Species  Continuity 
/da> 


EVi 


(10) 


(11) 


(12) 
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State 


I  a 


P2  =  P  2 


RT 


l2 

"2 


(13) 


2.3  SOLUTION  TECHNIQUE 


In  the  error  minimization  approach,  a  positive -definite  error  function 
is  defined  at  the  point  2;  vis:  . 


g2  =  2  [^1  EI  +  WII  EII  +  W 


III  EIII  +  WIV  EIV  + 


2>vi  EVi] 
i=l  J 


(14) 


where  the  Wj,  Wjj  W^TT,  W^,  WyT  are  scale  factors  qsed  to  ndndimensionalize 
the  corresponding  equations.  The  weighting  factors  were  based  on  the  ihtegral 
values  over  the  entire  n  station  and  applied  to  the  (n41)st  station 


max 


0  *  /  P^max- 

y=0 


f  max 

q  =  /  q  dy/y 

y=o 


max 


1 ' 


and 


1  1  1 

W  -  i W  -  -  W  W  -  - — 

I  -2-2’  II  "  -2-4  “  WIII»  WIV  — 2t6 


p  q 


p  q 


n-o-  Wvi  =  Wl 
P  q  1 


Now  the  point  2;  is  in  general,  any  interior  point  along  the  n+1  data  line  so  that 
a  positive -definite  error  function  for  the  entire  line  can  be  found  by  summing 
all  the  local  gs. 


- 1 

m=2 


n+1 


(15) 


The  object  is  to  reduce  G  to  an  acceptable  level  (zero  is  exact  but  impractical). 
This  reduction  may  be  made  by  the  method  of  steepest  descent.  Consider  that 
the  function  G  is  dependent  on  a  set  of  coordinates  £j;  then  from  the  calculus 
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/, 


where 


and 


dG  =  VG  • 


__  _  fV/1^  ^ 

T7r+  OVJ  y  ,  Caj  TT  *  Ou  ~ 

VG  ■  V'" wj*t. 


di  =  d|j  e6j  +d|2ee2  +  ...  +d6jeej 

I 

Now  fbr  a  given  step  length  |d|  |,  the  largest  change  or  payoff  in  G  occurs 
when  7G  and  Id  £  I  a^e  colinear.  Then  let 


so  that 


A?  vg 

-  |vojd$ 


VG  .  VG 

dG  -  |VG|  d* 


or 


dG  =  |VG|  d£ 

The  desired  change  in  G  is  from  its  current  level  to  zero,  i.e„ 

1 

dG  S  AG  =  0  -  Gk 


so  that 


and 


l||  =  A£  = 


-G' 


lvGkl 


aT  =  _Q]lZGk 
5  lvGkl2 


! 
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The  steepest  descent  recursion  formula  then  becomes 


|k+l 


Gk  VGk  .k 

k  ° 

VG  •  VG 


(17) 


where  a  step  modifier  6  has  been  included  to  control  the  descent  process. 

v 

An  unconditionally  stable  solution  can  be  obtained  by  altering  6  in  an  accept¬ 
able  fashion  if  the  merit  or  error  function  increases.  That  is  to  say  that  if 

G^+l  >  Gk, 

k  _ k+1 

then  the  step  modifier  6  is  halved  and  the  solution  tried  again  (£  is 

re-evaluated). 


To  proceed,  then  the  derivatives  of  G  with  respect  to  the  independent 
variables  (VG)  must  be  determined.  The  parameters  or  flow  properties  at 
station  n  are  fixed  so  that  G  must  be  differentiated  with  respect  to  (p,q,  9, 

T,a.)  for  all  points  on  the  n+1  data  line.  Take,  for  instance,  the  parameter 

*  t  s  t 

p  at  station  m'  on  the  n+1  data  line 


aa*£  _  a 

8G  _  \  '  fy?m.  n+1 
^m',  n+1  rn=2 

Because  of  the  three -point  influence  of  the  variable  pm  (the  station  subscript 
n+1  is  assumed  to  be  understood  from  this  point  on)  the  above  expression 
reduces  to 


80  ^m'-l  ,  88m',  88m'il  ... 

dp  /  dp  /  dp  r  dp  s 
Hm  Hm 

At  this  juncture  it  is  expedient  to  discuss  the  evaluation  of  the  derivatives  of 
gm,  with  respect  to  the  surrounding  points  and  reconstruct  the  above  gradient 
expression  later.  From  Eq.  (14)  it  can  be  seen  that  for  any  § 
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9g  8E.  9E_.  dE... 

-§f  ■  wiei1T4  wiiEn  ~W  *  wmEm  -# 


0E, 


IV 

+  Ww  EIV  -7j|-  +  Wv.  £,y.  — ij| 

i=l 


V'  8EVi 
i  2J  Evt 


(19) 


The  derivatives  of  E^  with  respect  to  the  variable  at  points  1,2  and  3 
are;  (all  derivatives  not  listed  are  zero) 


9EI  /00\  /8q\  q2sin02  .  q2 

apT  =  q2  W,  +  ifs).  +°  — 77“  +  AS 


9Ej 

W 


n  122.)  .  ^2  P2  sin92  . 

“  P2  *8n\  +AS  +  °  y,  +  Us' 


9Ej 

aeT 


p2q2 

2Ay 


9Ej  9Ej  P2q2  cosq; 

*V  =  °  ~Z  ' 


The  derivatives  of  Ejj  are: 


8E 


II 


dp. 


=  q 


aa  .  dIn 

2  8S  _  ’  9q, 


COS0- 


i  =  'liW'a 


aill  aa  +W2+jul.  !fiL  .  ,J_S_  +acJlh.\ 

an  as  +  s  +  ?>  p|  2  °2y2Ayj 


^2  "  K2  9S  2  ’  S  '  Ay2’  aq3 
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0ETT  sin0-  .  -  v 

11  ■  fe 


00, 


y2  '0n<2 


9EII  J_ 
dp~  AS 


The  derivatives  of  Ejjj  are 


0E 


III 


Op. 


=  q- 


(d e\ 

'OS'  ’  Oq 


III 

2 


OE 


III 


***»%  ~l2 


p_A 

AS 


OE 


III 

Opj 


1 

2Ay 


OE 


III 

0p3 


The  derivatives  of  E^^  are 


OE 


9p 


IV  _  r  /9T\  . 
r»-  '  q2  Cp  '  OS ’ 


OE 


IV 


Oq, 


Pz  cp  (H)2  "(Is  )2 


8EIV  jj_  (0 m  8EIV  8EIV  .  QM  CP  Sin92  (3T\ 

Oqj  “  Ay  '0n/2  0q3  ’  O02  y2  'On'. 


0ET„  q,  0E„,  P,q,  2  M  C 


TV  _  IV  _  2^2  p 

L-  A  c>  !vr  -  AS  + 


0p2  AS’  0T2 


Ay' 


OE 


_  a  r  \_J_  g  cose  I  iL  Y  c 

OTj  “  MCP)Ay2  2y2Ay(+Ay  Pf  '  9n  '2 


IV 
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j_i_  JCO,e2l  3Evi 

4Ly2'2>'2^r  \ 


+M. 

AS  A  2 


9EV1  |  i  a  c°.e2 1  9Eyi  y.  a  »n9z  /8a. . 

^ij  ‘  ‘"Ia,2  2  *2^  C  882  "  *2  '8n4 

2  2 

where  in  all  of  the  above  2Ay  =  -  y^,  while  Ay  =  (y^  -y^)  /4.  Also  in  the 

above  expressions,  the  derivatives  with  respect  to  p  have  been  taken  inde¬ 
pendently.  The  equations  of  state,  however,  relates  p  top,  T,a.,  The  chain 
rule  may  be  used  to  account  for  this  relationship.  For  instance 

dEIV 
9p2  9Pj 

p2  =  constant 
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where 


In  this  fashion,  the  gradient  vector  can  be  determined.  Components 
of  this  vector  may  exist  at  the  periphery  which  do  not  honor  or  satisfy  the 
boundary  conditions, 

2.4  BOUNDARY  CONDITIONS  ENFORCEMENT 

The  above  discussion  is  concerned  with  the  governing  equations  of 
motion  controlling  the  interior  region  of  the  flow  field.  In  addition  to  these 
relationships,  information  must  be  supplied  at  the  boundaries.  These  con¬ 
ditions  select  from  a  large  number  of  possible  solutions,  that  solution  which 
is  appropriate  to  the  particular  problem.  These  conditions  are  not  always 
easy  to  specify  physically  or  mathematically.  Fortunately,  in  the  parabolic 
mixing  analysis  the  boundary  conditions  are  straightforward. 

At  large  lateral  distances,  the  flow  variables  must  return  to  the  free- 
stream  values.  If  the  data  points  were  initialized  to  these  values  then  one 
need  only  require  that 


no  no  no  no  no 

<AJ  _  (A-J  CAJ  CAJ  _  CAJ  ~ 

dp  ~  9q  =  90  =  9T  "  acr  = 

At  the  inner  wall,  however,  the  situation  is  somewhat  more  complicated. 
Although  the  equations  (for  the  axisymmetric  case)  have  been  written  with 
implied  symmetry,  there  is  really  nothing  in  the  numerical  analog  to 
guarantee  this.  Because  of  symmetry  certain  functions  are  known  to  be 
symmetric  about  the  centerline  while  others  are  antisymmetric.  For 
instance  the  axial  velocity  component,  density,  temperature,  and  species 
are  of  equal  value  above  and  below  the  centerline  while  the  flow  angle  is 
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equal  in  magnitude  but  opposite  in  sign.  Therefore,  for  all  but  the  flow 
angle  (£  stands  for  any  variable)  the  boundary  condition  is 

91/ 9y  =  0 

A  second -order  Taylor  series  expansion  in  this  region  in  which  point  1  is  on 
the  axis,  2  is  just  above  the  axis  and  3  is  the  point  above  2.  So  that 

6  = it  .it 

61  3  ’2  3  s3 


In  taking  derivatives  previously,  point  1  was  considered  independently 
of  points  2  and  3.  The  symmetry  condition,  however,  causes  a  dependency. 
Again  using  the  chain  rule  this  effect  may  be  accounted  for  by  writing 


and 


3G 


*1 

ai2 


8G 


dC 

91, 


*1 

*3 


When  the  above  operation  has  been  completed,  the  gradient  vector  has  been 
correctly  determined,  and  the  integration  using  the  recursion  formula,  Eq. 
(17),  can  be  performed. 
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Section  3 
RESULTS 

A  typical  turbojet  engine/atmospheric  mixing  case  is  shown  as  a  sample 
of  the  program  predictions.  The  case  analyzed  is  a  balanced  jet  case  in  which 
the  jet  velocity  is  1800  ft/sec  while  the  freestream  velocity  is  880  ft/sec.  The 
jet  temperature  is  900°K  and  the  freestream  is  300°K.  A  constant  eddy  vis¬ 
cosity  model  was  used  with  the  eddy  viscosity  chosen  to  be  0.01.  For  this 
example  case  both  streams  were  assumed  to  have  a  molecular  weight  of  28.02 
and  frozen  (one  species)  flow  was  assumed. 

Figure  1  shows  the  velocity  distribution  versus  lateral  distance  across 
the  mixing  zone.  Note  that  the  inward  mixing  rate  is  much  higher  than  the 
outer  rate.  The  velocity  profiles  seem  to  be  affected  somewhat  by  the  duct 
wall  or  outer  boundary  condition,  and  this  may  affect  the  inward  rate.  In  a 
subsequent  figure,  however,  it  will  be  shown  that  the  temperature  spread  does 
not  act  in  the  same  fashion. 

The  velocity  distributions  are  shown  for  X  =  0,  1, 2,  4,  8  ft  from  the  jet 
exit  plane.  Figure  2  illustrates  the  axial  decay  of  the  jet  centerline  velocity. 
The  initial  decay  rate  appears  to  be  much  too  fast  (eddy  viscosity  is  too  large) 
causing  the  potential  core  to  disappear  almost  immediately  while  the  decay 
rate  farther  downstream  is  more  realistic  indicating  that  the  eddy  viscosity 
chosen  is  appropriate  at  large  distances. 

Figure  3  presents  the  temperature  distribution.  The  inner  and  outer 
spread  rates  seem  to  be  nearly  equal.  This  is  not  the  same  behavior  shown 
by  the  velocity  distribution  and  would  indicate  that  the  duct  wall  boundary 
conditions  are  not  affecting  the  analysis  as  previously  conjectured.  There 
is  no  valid  reason  to  assume  that  both  spread  rates  should  be  the  same  nor 
should  the  temperature  mixing  rate  be  the  same  as  that  of  the  velocity. 
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200 


A  X  =  0  (ft) 
O  X  =  1 
□  X  =  2 
Q  X  =  4 
O  X  =  8 


Velocity  Distribution  vs  Lateral  Distance  for  Typical  Turbojet 
Mixing  Case 


Axial  Distance  (ft) 


Fig.  2  -  Jet  Centerline  Velocity  Decay 
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Temperature  (°K) 


Fig.  3  -  Temperature  Distribution  vs  Lateral  Distance  for  Typical 
Turbojet  Mixing 


LOCKHEED  -  HUNTSVILLE  RESEARCH  &  ENGINEERING  CENTER 


LMSC-HREC  D306224 


The  analysis,  being  implicit,  does  not  suffer  from  a  stability  problem 
so  that  large  axial  steps  may  be  taken  and  profile  definitions  may  be  obtained 
at  large  axial  distances  with  nominal  computational  expense. 
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Section  4 
CONCLUSIONS 

The  results  obtained  indicate  that  the  computer  program  will  produce 
meaningful  results  in  a  small  amount  of  computer  time.  More  analysis  is 
required,  however,  to  assess  the  validity  of  the  calculations  for  imbalanced 
jets  and  reacting  mixtures. 

In  conclusion  it  is  pointed  out  that  repetitive  temperature  peaks  are 
within  the  basic  capability  of  the  computer  program  for  all  speed  regimes. 
The  formation  of  Mach  disks,  however,  being  elliptic  in  nature  will  not  be 
adequately  predicted  although  it  is  possible  that  one  of  the  several  approxi¬ 
mate  Mach  disk  location  theories  may  be  used  to  form  the  disk  and  that 
the  calculation  could  then  be  resumed  in  the  parabolic  fashion. 
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Appendix 

COAXIAL  JET  MIXING  WITH  LATERAL  PRESSURE 
GRADIENTS  COMPUTER  PROGRAM 


Appendix 

1.0  INTRODUCTION 

The  computer  program  for  "Coaxial  Jet  Mixing  with  Lateral  Pressure 
Gradients"  is  written  in  FORTRAN  V  for  the  CDC  6600  computer  at  the  Army 
Missile  Command  Computation  Center  at  Redstone  Arsenal.  Little  or  no  dif¬ 
ficulty  is  anticipated  in  program  utilization  at  other  facilities  with  different 
machines.  This,  appendix  is  concerned  with  the  computer-oriented  aspects  of 
the  anal/sis  as  opposed  to  the  theoretical  considerations  which  are  discussed 
n  the  main  body  of  the  report. 

2.0  DISCUSSION  ■ 

2.1  Subroutine  Description 

i 

The  program  is  arranged  in  logical  subdivisions  or  subroutines.  The 
routines,  other  than  system  routines,  are-  MAIN,  OUTPUT,  LOGIC,  EMT, 
and  INOUT.  The  primary  function  of  each  routine  is  given  below: 

I 

MAIN  —  In  this  routine  the  basic  data  and  run  control  information  are 
read  and  the  initial  data  surface  for  the  forward  marching  solution  is  prepared. 
The  routine  also  computes  thermodynamic  properties  and  species  production 
rate.  'An  initial  guess  for  the  solution  at  the  next  (downstream)  surface  is  pre¬ 
pared  and  the  flow  solution  subroutines  are  called.  Upon  termination  of  the 
flow  solution  for  the  new  surface,  control  is  returned  to  MAIN. 

I 

1  » 

INOUT  —  This  subroutine  outputs  the  input  values 

) 

LOGIC  —  This  subroutine  controls  the  descent  process.  Based  on  the 
gradient  values  determined  in  EMT  a  step  length  of  each  variable  in  the  field 
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is  determined  and  the  next  approximation  to  the  flow  variables  is  made.  In 
the  event  that  the  merit  function  increases  rather  than  decreases  the  step 
length  is  halved.  A  cyclic  process  is  used  rather  than  steepest  descent  in 
that  only  velocity  changes  are  considered  until  a  successful  step  is  made, 
then  temperature,  flow  angle,  density,  and  species  are  processed.  After  a 
cycle  has  been  completed  it  is  repeated  until  the  maximum  number  of  itera¬ 
tions  have  been  performed.  Control  is  returned  to  MAIN  and  the  next  march¬ 
ing  step  is  taken. 

EMT  —  The  merit  function  and  its  derivatives  are  calculated  in  this 
subroutine.  See  the  main  body  of  the  report  for  a  discussion  of  these 
calculations. 

OUTPUT  —  This  routine  outputs  the  current  axial  step  flow  properties. 
The  overall  program  flow  is  given  below: 


MAIN  INOUT 


A-2 


LOCKHEED  -  HUNTSVILLE  RESEARCH  &  ENGINEERING  CENTER 


LMSC-HREC  D306224 


2.2  INPUT  INSTRUCTIONS 


CARD  1 


CARD  2 


CARD  3 


FORMAT  (12A6) 

FIELDS  1-12  Any  Hollerith  information,  heading 
printed  atop  each  page 


FORMAT 
FIELD  1 
FIELD  2 

FIELD  3 

FIELD  4 


FIELD  5 


FIELD  6 

FIELD  7 

FIELD  8 
FIELD  9 
FIELD  10 


(1615)  -  Run  control  parameters 

Number  of  radial  stations  in  jet  stream 

Number  of  radial  stations  in  atmospheric 
stream 

Number  of  input  cards  with  jet  properties 
equal  1  for  uniform  jet  setup  equal  FIELD  1 
if  distribution  across  jet  desired 

Number  of  input  cards  with  atmospheric 
properties  equal  1  for  uniform  atmos¬ 
pheric  stream  equal  FIELD2  if  distri¬ 
bution  desired 

Integer  selecting  type  of  eddy  viscosity 
model 

1  -  laminar  (Sutherland's)  model 

2  -  constant  viscosity  model 

3  -  Ferri  viscosity  model 

4  -  Ting-Libby  viscosity  model 

5  -  Ting-Libby  after  mixing  region 

intersects  axis 

6  -  GASL  mixing  width  model 

Number  of  species  present  across  entire 
region 

0  two-dimensional,  1  for  axisymmetric 
analysis 

Number  of  relaxation  steps/print  step 

Maximum  number  of  iterations 

Number  of  reactions  (zero  if  frozen 
calculation) 


FIELD  11-16 
FORMAT  7E10 
FIELD  1 
FIELD  2 
FIELD  3 
FIELD  4 


Not  presently  used 

.6 

Axial  coordinate  of  initial  surface 
Axial  step  length 
Maximum  axial  distance 
Initial  jet  radius 
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CARD  3 


CARD 

GROUP(S)  4 

CARD  1 


CARD(S)  2 
CARD 

GROUP(S)  5 

CARD  1 


CARDS(2) 

CARD 

GROUP(S)6 


FORMAT 
FIELD  'j 
FIELD  6 
FIELD  7 


7E10.6 

Initial  outer  edge  of  atmospheric  stream 
Normalizing  dimension  (usually  jet  radius) 
Eddy  viscosity  or  coefficient  as  appropriate 


FORMAT  (7E10.6)  -  As  many  card  groups  as  in  FIELD3  of 
CARD  2  -  the  Lrst  card  group  pertains  to 
jet  centerline 


FIELD  1  Pressure  (atm)  of  jet 

FIELD  2  Temperature  (°K)  of  jet 

FIELD  3  Velocity  (ft/sec)  of  jet 

FIELD  4  Angle  (deg)  of  jet 

FIELD  5  Radial  coordinate  associated  with  these 

properties 

FIELDS(l-NS)  Species  mole  fractions  at  this  location 


FORMAT  (7E 

FIELD  1 
FIELD  2 
FIELD  3 
FIELD  4 
FIELD  5 

FIELD(l-NS) 


10.6)  -  As  many  cards  as  in  FIELD4  of 
CARD  2  -  the  first  card  pertains  to  the 
freestream  edge  just  above  the  jet 

Pressure  (atm)  of  atmosphere 

Temperature  (°R)  of  atmosphere 

Velocity  (ft/sec)  of  atmosphere 

Angle  (deg)  of  atmosphere 

Radial  coordinate  associated  with  these 
properties 

Species  mole  fractions  at  this  location 


FORMAT  (A6/4E16.8/4E16.8/4E16.8)  -  Species  infor¬ 
mation.  1  card  group  for  each  species 


FIELD  1 
FIELD  2 
FIELD  3 
FIELD  4 
FIELD  5 
FIELD  6-13 


Species  name 
Molecule  weight 
C 

g 


A-4 


LOCKHEED  HUNTSVILLE  RESEARCH  &  ENGINEERING  CENTER 


LMSC-HREC  D306224 


where 


and 


Cg  =  reference  entropy,  =  reference  enthalpy 
C  =  }  C.  Tl 

P  l 

i=-4 


C  j  —  C  ^2  ~  ^  3’ 


CARD(S)  7  FORMAT 

FIELD  1 
2 

3 

4 

5 

6 

7 

8 

9 

10 


(A6,  IX,  A6,  8X,  A6,  IX,  A6,1X,A6,  7X,  12,11, 
E8.2,  F4.1,  F9.1)  read  in  NR  such  cards  in 
any  order 

Name  of  species  A  A6 

Name  of  species  B  A6 

Name  of  species  C  A6 

Name  of  species  D  (or  M)  A6 

Name  of  species  E  (or  M)  A6 

Reaction  type  1-10  (Ref.  1)  12 

Rate  constant  type  (Ref.  1)  II 

Pre-exponentional  factor  (cm- 
particle-sec)  E8.2 

Temperature  exponent  F4.1 

Activation  energy  cal/mole  F9.1 


The  above  input  defines  the  reaction  and  its  rate.  The  reaction  equa 
tion  is 

A+B  — ►C+D+E 

while  the  rate  equation  is  k  =  aT"*5  i~c/^  . 
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SAMPLE  OUTPUT 
Page  3 

(Note:  Only  one  species  present) 


x*  o. 


PT  Y/* 


1  0.  000  00 
2  .  10000 

3 

4  .  30000 

5  .40000 

6  .50000 

7  .60000 

8  .  70000 

9  .  48098- 

10  .90000 

11  1.  00000 
12  1.  10000 

13  1.20000 

14  1.  30.0  00 

15  1 *-400-80 — 

16  1.50030 

17  1.60000 

18  1.700  30 

19  1.80000 

20  1.90000 

21  2.00000 


FEET 


N2 


1.08  36UE-+86 
1.0300 GE+Oo 

1.00 J0GE  +  00 
1.00-08  0E  +  Q6 
1.0009  OE ♦  0 C 
l.OOOOOE+OO 
1.3000  GE  +  00 

1 • 0903  GE  +  OC 
l.OOOOGE+Ot* 
l.OOCOOE+OC 
1.0000  GE+09 
1.3000  OE  *-00 
l.OOOOOE+OO 
1 . 0303  0E+30 
1. OOOOGE+Ofl 
1.3000 CE+QC 
1 .00088E+8G 
l.OOOOOE+OO 
-i .  o4oooo+oo 
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